In [93]:
from astropy.io import fits
from astropy.table import Table
import numpy as np
import matplotlib.pyplot as plt
from astropy.stats import sigma_clipped_stats
from photutils import datasets
from photutils import CircularAperture
from photutils import DAOStarFinder

hdu = datasets.load_star_image('quadRU.fits')    
data = hdu.data#[0:401, 0:401]    
mean, median, std = sigma_clipped_stats(data, sigma=3.0)    
print((mean, median, std)) 
(3770.4896097401975, 3744.0, 237.21215226532127)
In [94]:
#Plotting the image
fig, ax = plt.subplots(1,2)

fig.set_size_inches(16,12)

ax[0].imshow(data,cmap=plt.cm.gray) #change colour to gray => cmap=plt.cm.gray

ax[1].imshow(data[401:801, 401:801] ,cmap=plt.cm.gray)

print("The number of pixels is {0}".format(data.size))
print("The mean is {0} and the STD ".format(np.mean(data)),format(np.std(data)))
The number of pixels is 1121481
The mean is 3988.5913154123878 and the STD  1069.1249436838032
In [95]:
#STEP 6
daofind = DAOStarFinder(fwhm=3.0, threshold=27.7*std)    #Full width at half maximum(fwhm) // good value = 27.7
sources = daofind(data)# - median)    
for col in sources.colnames:    
    sources[col].info.format = '%.8g'  # for consistent table output
sources.sort('xcentroid')
sources
Out[95]:
Table length=208
idxcentroidycentroidsharpnessroundness1roundness2npixskypeakfluxmag
int64float64float64float64float64float64int64float64float64float64float64
1700.91894685840.895150.300171390.423740310.73439404250126451.4511468-0.40427836
1321.3110953603.525030.3070012-0.175347050.8687376250125861.0172121-0.018528787
481.3819622194.148960.34887691-0.0403594720.91937352250124391.0000395-4.2841351e-05
19413.709729966.332180.37653395-0.061260617-0.0047948456250125511.0262129-0.028093697
2814.85393104.588580.41325279-0.157549590.0014790695250121761.078323-0.081872155
14315.308049662.256540.404104940.092247523-0.21561658250127341.0472938-0.050171342
20126.026207990.829510.40973918-0.12748799-0.1263963250122891.1639493-0.16483515
4131.916545171.798050.409471790.073722424-0.097077158250126961.1165765-0.11972121
6335.353299267.108790.40362557-0.1011558-0.10080115250124981.0806916-0.084254463
2137.11509983.2903520.553584020.0189599820.032717845250115301.007856-0.00849618
.................................
151974.21449695.409980.48341372-0.313254910.011266162250119351.0067901-0.0073473202
5981.8168232.8097470.49815515-0.038007005-0.15331728250116821.0966799-0.10019975
1111014.2458456.118190.382510060.0754759580.027730015250125291.0789151-0.082468222
1251023.7488564.62620.440421570.30361012-0.042326908250116001.0066716-0.007219594
1141039.4274492.844680.438376810.031428091-0.24080226250123201.0214756-0.023070039
271041.6488101.708370.454023990.19568416-0.058768625250118251.0237259-0.025459175
351043.7634119.084710.494115840.031467285-0.10235053250121431.1573171-0.15863093
1381044.666632.901660.53714543-0.1248506-0.12651939250118951.0521495-0.055193654
1671044.7256804.824160.53046292-0.11265257-0.044356247250115941.0156725-0.016884211
961048.2053380.564770.41206942-0.19084793-0.080891852250124451.0497809-0.052746648
In [ ]:
 
In [96]:
import matplotlib.pyplot as plt
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from photutils import CircularAperture

plt.figure(figsize=(20,20))
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=4.)

norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(data, cmap='Greys', origin='lower', norm=norm)

apertures.plot(color='blue', lw=1.5, alpha=0.5)
In [ ]:
 
In [ ]:
 
In [97]:
#STEP 8
import matplotlib.pyplot as plt
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from photutils import CircularAperture
from photutils import aperture_photometry

plt.figure(figsize=(20,20))
positions = (sources['xcentroid'], sources['ycentroid'])

apertures = CircularAperture(positions, r=4.)
phot_table = aperture_photometry(data, apertures, method='exact',subpixels=5)

norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
apertures.plot(color='blue', lw=1.5, alpha=0.5)
In [98]:
phot_table = aperture_photometry(data, apertures, method='exact',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g'  # for consistent table output
print(phot_table)
 id      xcenter            ycenter       aperture_sum
           pix                pix                     
--- ------------------ ------------------ ------------
  1 0.9189468499910883  840.8951457248296    235087.24
  2  1.311095259132758  603.5250306712844    387140.03
  3 1.3819621776452984  194.1489634549567    372487.51
  4 13.709729099697494  966.3321753368325    298805.15
  5 14.853929564441659  104.5885806490168    275656.75
  6 15.308049034659224  662.2565373858495    320886.39
  7 26.026207222819796  990.8295092551749    271741.42
  8  31.91654474923815 171.79804519605577    283941.14
  9 35.353299330115014  267.1087864002644    285376.66
 10  37.11509862938315  83.29035164102501    249927.49
...                ...                ...          ...
199  974.2144947330132  695.4099754883719    265665.99
200  981.8168232536541  32.80974670705418    247157.03
201   1014.24581863924  456.1181928875807    293596.48
202 1023.7488485932108  564.6262045991909    266480.28
203 1039.4273600411575 492.84467803095777    287349.87
204  1041.648750654047  101.7083747902434    267495.91
205  1043.763408945062 119.08471405894595    254575.88
206  1044.666008741653  632.9016640068996    256337.64
207  1044.725609535972  804.8241634613814    253547.55
208 1048.2053425810136 380.56477025247375    292795.78
Length = 208 rows
In [99]:
phot_table = aperture_photometry(data, apertures, method='subpixel',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g'  # for consistent table output
print(phot_table)
 id      xcenter            ycenter       aperture_sum
           pix                pix                     
--- ------------------ ------------------ ------------
  1 0.9189468499910883  840.8951457248296     235553.8
  2  1.311095259132758  603.5250306712844    388696.92
  3 1.3819621776452984  194.1489634549567    372579.88
  4 13.709729099697494  966.3321753368325    299362.88
  5 14.853929564441659  104.5885806490168    275714.36
  6 15.308049034659224  662.2565373858495    321257.12
  7 26.026207222819796  990.8295092551749     271478.6
  8  31.91654474923815 171.79804519605577    283271.96
  9 35.353299330115014  267.1087864002644    285431.44
 10  37.11509862938315  83.29035164102501    250929.24
...                ...                ...          ...
199  974.2144947330132  695.4099754883719    265225.84
200  981.8168232536541  32.80974670705418    246734.24
201   1014.24581863924  456.1181928875807    293849.52
202 1023.7488485932108  564.6262045991909    266887.52
203 1039.4273600411575 492.84467803095777    287735.48
204  1041.648750654047  101.7083747902434    267436.12
205  1043.763408945062 119.08471405894595    254792.08
206  1044.666008741653  632.9016640068996    256692.12
207  1044.725609535972  804.8241634613814    253284.56
208 1048.2053425810136 380.56477025247375    292364.44
Length = 208 rows
In [100]:
phot_table = aperture_photometry(data, apertures, method='center',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g'  # for consistent table output
print(phot_table)
 id      xcenter            ycenter       aperture_sum
           pix                pix                     
--- ------------------ ------------------ ------------
  1 0.9189468499910883  840.8951457248296       226283
  2  1.311095259132758  603.5250306712844       404626
  3 1.3819621776452984  194.1489634549567       371062
  4 13.709729099697494  966.3321753368325       298020
  5 14.853929564441659  104.5885806490168       274731
  6 15.308049034659224  662.2565373858495       326327
  7 26.026207222819796  990.8295092551749       266936
  8  31.91654474923815 171.79804519605577       283197
  9 35.353299330115014  267.1087864002644       288381
 10  37.11509862938315  83.29035164102501       252774
...                ...                ...          ...
199  974.2144947330132  695.4099754883719       260311
200  981.8168232536541  32.80974670705418       254090
201   1014.24581863924  456.1181928875807       296865
202 1023.7488485932108  564.6262045991909       265634
203 1039.4273600411575 492.84467803095777       286708
204  1041.648750654047  101.7083747902434       266108
205  1043.763408945062 119.08471405894595       253822
206  1044.666008741653  632.9016640068996       259520
207  1044.725609535972  804.8241634613814       256679
208 1048.2053425810136 380.56477025247375       291884
Length = 208 rows
In [137]:
#STEP 10 - 11
from photutils import CircularAnnulus

apertures4 = CircularAperture(positions, r=3.)
annulus_apertures4 = CircularAnnulus(positions ,r_in=6 ,r_out=9)
apers4 =  [apertures4, annulus_apertures4]

phot_table4 = aperture_photometry(data, apers4)

for col in phot_table4.colnames:
    phot_table4[col].info.format = '%.8g'  # for consistent table out
print(phot_table4)

#The aperture_sum_0 column refers to the first aperture (circular)
#the aperture_sum_1 column refers to the second aperture (annulus)
 id  xcenter    ycenter  aperture_sum_0 aperture_sum_1
       pix        pix                                 
--- ---------- --------- -------------- --------------
  1 0.91894685 840.89515      175857.43       299406.7
  2  1.3110953 603.52503       283355.7      342294.06
  3  1.3819622 194.14896      280215.01      317289.96
  4  13.709729 966.33218      206998.37      524160.52
  5   14.85393 104.58858      187510.09      533347.69
  6  15.308049 662.25654      210624.24      594384.47
  7  26.026207 990.82951      183965.79      521889.02
  8  31.916545 171.79805      194901.82       518105.3
  9  35.353299 267.10879      195810.43      514602.35
 10  37.115099 83.290352         160962      535694.26
...        ...       ...            ...            ...
199  974.21449 695.40998      178183.94      537045.49
200  981.81682 32.809747      162077.25      522570.44
201  1014.2458 456.11819      200126.95      548544.89
202  1023.7488  564.6262      177732.56      544490.79
203  1039.4274 492.84468      194633.15      537354.49
204  1041.6488 101.70837      178287.04      524084.66
205  1043.7634 119.08471      168338.84      520448.17
206   1044.666 632.90166      166209.24      550289.95
207  1044.7256 804.82416      164841.21      536903.28
208  1048.2053 380.56477      198785.18      547217.71
Length = 208 rows
In [138]:
bkg_mean = phot_table4['aperture_sum_1'] / annulus_apertures.area()
In [139]:
bkg_sum = bkg_mean * apertures4.area()

final_sum = phot_table4['aperture_sum_0'] - bkg_sum
phot_table4['residual_aperture_sum'] = final_sum
phot_table4['residual_aperture_sum'].info.format = '%.8g'  # for consistent table output
print(phot_table4['residual_aperture_sum'])  
residual_aperture_sum
---------------------
            154300.15
            258710.53
            257370.14
            169258.81
            149109.06
            167828.56
            146389.78
            157598.24
            158759.06
            122392.01
                  ...
            120263.39
            139516.67
            124452.18
            160631.72
            138529.23
            155943.63
            140552.95
            130866.57
            126588.36
            126184.18
            159385.51
Length = 208 rows
In [140]:
from photutils import CircularAperture, CircularAnnulus
from astropy.visualization import simple_norm

norm = simple_norm(data, 'sqrt', percent=99)
plt.figure(figsize=(30,30))
plt.imshow(data, norm=norm)
apertures4.plot(color='white', lw=2)
annulus_apertures4.plot(color='red', lw=2)
In [144]:
import matplotlib.pyplot as plt

plt.figure(figsize=(15,15))
plt.scatter(phot_table4['residual_aperture_sum'], phot_table4['aperture_sum_0'])
plt.xlabel("Circular Apertures")
plt.ylabel("Background")
plt.xlim(100000,200000)
plt.ylim(150000,240000)
Out[144]:
(150000, 240000)
In [150]:
pos = Table(names=['x_0', 'y_0'], data=[sources['xcentroid'],
                                        sources['ycentroid']])
In [180]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from photutils.datasets import (make_random_gaussians_table,
                                make_noise_image,
                                make_gaussian_sources_image)
from photutils.psf import (IterativelySubtractedPSFPhotometry,
                           BasicPSFPhotometry)
from photutils import MMMBackground
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.detection import DAOStarFinder
from photutils.detection import IRAFStarFinder
from astropy.table import Table
from astropy.modeling.fitting import LevMarLSQFitter
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.background import MMMBackground, MADStdBackgroundRMS
from photutils.psf import IterativelySubtractedPSFPhotometry
from astropy.stats import gaussian_sigma_to_fwhm
sigma_psf = 2.0
bkgrms = MADStdBackgroundRMS()
std = bkgrms(data)
iraffind = IRAFStarFinder(threshold=3.5*std,
                          fwhm=sigma_psf*gaussian_sigma_to_fwhm,
                          minsep_fwhm=0.01, roundhi=5.0, roundlo=-5.0,
                          sharplo=0.0, sharphi=2.0)

daogroup = DAOGroup(2.0*sigma_psf*gaussian_sigma_to_fwhm)
mmm_bkg = MMMBackground()
fitter = LevMarLSQFitter()
psf_model = IntegratedGaussianPRF(sigma=sigma_psf)



photometry = BasicPSFPhotometry(group_maker=daogroup,
                                bkg_estimator=mmm_bkg,
                                psf_model=psf_model,
                                fitter=LevMarLSQFitter(),
                                fitshape=(11,11))
result_tab = photometry(image=data, init_guesses=pos)
result_tab
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Out[180]:
Table length=208
x_0y_0flux_0idgroup_idx_fity_fitflux_fitflux_uncx_0_uncy_0_unc
float64float64float64int64int64float64float64float64float64float64float64
0.91894685840.89515104622.9867255080311-0.3173792241660808840.8844099643336228890.4305793325811388.2315068218160.149680729652099280.09072328854096326
1.3110953603.52503261117.40290235524221.300569260807128603.1321144660371320870.12180474089054.566572061230.093834157882552920.0787160379628428
1.3819622194.14896237618.21822195407331.3695821448723215194.1425495250382303854.0648381155582.7237398865950.060730666006209810.05136082396918347
13.709729966.33218117309.864778372924413.712264141803514966.2986230219999159191.537017787226206.9662153025380.111510329058488340.11157241528485883
14.85393104.5885893019.761385789025514.891850631606813104.72455897028448131748.695882852046487.7747911223060.140840375690643060.14086124526028662
15.308049662.25654162433.684013514756616.56388179248046661.5843508658699184976.2205967409514269.970784404530.222977976701880380.220723213395023
26.026207990.8295189210.758820951257726.02352000642823990.787169671028126253.197130868926245.7214753191970.141453052704833850.14152292535230906
31.916545171.79805100658.136800201578831.917581972682235171.8343283450299141613.772637973426448.4237872656230.130232376866731340.13023927840828783
35.353299267.10879103527.162757632839935.350323440435055267.06805584474233143578.564579712876366.092452122290.126844912519184420.12681320427197512
37.11509983.29035269293.6427486517101037.10641610557479583.215215885896192957.740369714485166.6423027881630.158964808325431640.15897906232912956
.................................
974.21449695.4099883716.48918003912199195974.1577370810207695.4605601376018119046.200686188526389.6349090896630.153539453159302150.15358879510160484
981.8168232.80974764075.48566281513200196981.814758405302732.8093077725045694197.524286495865928.2168695678880.180037126872748530.17996729945719703
1014.2458456.11819115649.236801013932011971014.2525036848132456.0948214803198151671.880558922546114.4576026766710.115347106897731450.1152792350735681
1023.7488564.626285052.004606317342021981023.7318537578201564.7185501263821117954.23027648245902.8155710030130.143162000886637440.14312117353528225
1039.4274492.84468107735.284469785632031991039.482610771519492.79794173036663143473.74462374936066.0738626865040.121038451046032240.12090424187403902
1041.6488101.7083786397.264132560682042001041.6762768184842101.69553365408699119421.107676376116064.9479538548780.14530710723261630.14524493521378323
1043.7634119.0847171455.857626797122052011043.7819088676713119.1137619317121103288.88710708716018.5552372328580.166699774403734490.16661564847624116
1044.666632.9016677439.811920777252062021044.7029909726716632.9280085526168102270.056992100845505.3668660392390.154002246550849260.153935963968129
1044.7256804.8241672539.80436439012072031044.728165492963804.855394634385798589.338260883425259.2658936537110.152598180860739760.1525560134208171
1048.2053380.56477114486.678459136972082041048.1557223219777380.56131251673327149602.081117475735893.0516125452790.11269841743514050.11268198008765348
In [185]:
#norm = ImageNormalize(strecth=SqrtStrecth())

plt.figure(figsize=(15,15))
plt.scatter(phot_table4['aperture_sum_0'], result_tab['flux_fit'])

plt.xlabel("Aperture Sum")
plt.ylabel("PSF Value")

plt.xlim(150000,240000)
#plt.ylim(150000,240000)
Out[185]:
(150000, 240000)
In [169]:
plt.figure(figsize=(30,30))
plt.subplot(1, 2, 1)

plt.imshow(data, cmap='viridis', aspect=1,
           interpolation='nearest', origin='lower')
plt.title('Simulated data')
plt.colorbar(orientation='horizontal', fraction=0.046, pad=0.04)
plt.subplot(1 ,2, 2)
plt.imshow(residual_image, cmap='viridis', aspect=1,
           interpolation='nearest', origin='lower')
plt.title('Residual Image')
plt.colorbar(orientation='horizontal', fraction=0.046, pad=0.04)
Out[169]:
<matplotlib.colorbar.Colorbar at 0xb1b691128>
In [ ]: